home *** CD-ROM | disk | FTP | other *** search
/ Amiga Format CD 41 / Amiga Format CD41 (1999-06)(Future Publishing)(GB)[!][issue 1999-07].iso / -seriously_amiga- / programming / other / scm / slib / randinex.scm < prev    next >
Text File  |  1999-04-19  |  5KB  |  126 lines
  1. ;;;"randinex.scm" Pseudo-Random inexact real numbers for scheme.
  2. ;;; Copyright (C) 1991, 1993, 1999 Aubrey Jaffer.
  3. ;
  4. ;Permission to copy this software, to redistribute it, and to use it
  5. ;for any purpose is granted, subject to the following restrictions and
  6. ;understandings.
  7. ;
  8. ;1.  Any copy made of this software must include this copyright notice
  9. ;in full.
  10. ;
  11. ;2.  I have made no warrantee or representation that the operation of
  12. ;this software will be error-free, and I am under no obligation to
  13. ;provide any services, by way of maintenance, update, or otherwise.
  14. ;
  15. ;3.  In conjunction with products arising from the use of this
  16. ;material, there shall be no use of my name in any advertising,
  17. ;promotional, or sales literature without prior written consent in
  18. ;each case.
  19.  
  20. ;This file is loaded by random.scm if inexact numbers are supported by
  21. ;the implementation.
  22.  
  23. ;;; Sphere and normal functions corrections from: Harald Hanche-Olsen
  24.  
  25. ;;; Generate an inexact real between 0 and 1.
  26. (define random:uniform1
  27.                     ; how many chunks fill an inexact?
  28.   (do ((random:chunks/float 0 (+ 1 random:chunks/float))
  29.        (smidgen 1.0 (/ smidgen 256.0)))
  30.       ((or (= 1.0 (+ 1 smidgen)) (= 4 random:chunks/float))
  31.        (lambda (state)
  32.      (do ((cnt random:chunks/float (+ -1 cnt))
  33.           (uni (/ (random:chunk state) 256.0)
  34.            (/ (+ uni (random:chunk state)) 256.0)))
  35.          ((= 1 cnt) uni))))))
  36.  
  37.  
  38. ;;@args
  39. ;;@args state
  40. ;;Returns an uniformly distributed inexact real random number in the
  41. ;;range between 0 and 1.
  42. (define (random:uniform . args)
  43.   (random:uniform1 (if (null? args) *random-state* (car args))))
  44.  
  45.  
  46. ;;@args
  47. ;;@args state
  48. ;;Returns an inexact real in an exponential distribution with mean 1.  For
  49. ;;an exponential distribution with mean @var{u} use
  50. ;;@w{@code{(* @var{u} (random:exp))}}.
  51. (define (random:exp . args)
  52.   (- (log (random:uniform1 (if (null? args) *random-state* (car args))))))
  53.  
  54.  
  55. ;;@args
  56. ;;@args state
  57. ;;Returns an inexact real in a normal distribution with mean 0 and
  58. ;;standard deviation 1.  For a normal distribution with mean @var{m} and
  59. ;;standard deviation @var{d} use
  60. ;;@w{@code{(+ @var{m} (* @var{d} (random:normal)))}}.
  61. (define (random:normal . args)
  62.   (let ((vect (make-vector 1)))
  63.     (apply random:normal-vector! vect args)
  64.     (vector-ref vect 0)))
  65.  
  66.  
  67. ;;; If x and y are independent standard normal variables, then with
  68. ;;; x=r*cos(t), y=r*sin(t), we find that t is uniformly distributed
  69. ;;; over [0,2*pi] and the cumulative distribution of r is
  70. ;;; 1-exp(-r^2/2).  This latter means that u=exp(-r^2/2) is uniformly
  71. ;;; distributed on [0,1], so r=sqrt(-2 log u) can be used to generate r.
  72.  
  73. ;;@args vect
  74. ;;@args vect state
  75. ;;Fills @1 with inexact real random numbers which are independent
  76. ;;and standard normally distributed (i.e., with mean 0 and variance 1).
  77. (define (random:normal-vector! vect . args)
  78.   (let ((state (if (null? args) *random-state* (car args)))
  79.     (sum2 0))
  80.     (let ((do! (lambda (k x)
  81.          (vector-set! vect k x)
  82.          (set! sum2 (+ sum2 (* x x))))))
  83.       (do ((n (- (vector-length vect) 1) (- n 2)))
  84.       ((negative? n) sum2)
  85.     (let ((t (* 6.28318530717958 (random:uniform1 state)))
  86.           (r (sqrt (* -2 (log (random:uniform1 state))))))
  87.       (do! n (* r (cos t)))
  88.       (if (positive? n) (do! (- n 1) (* r (sin t)))))))))
  89.  
  90.  
  91. ;;; For the uniform distibution on the hollow sphere, pick a normal
  92. ;;; family and scale.
  93.  
  94. ;;@args vect
  95. ;;@args vect state
  96. ;;Fills @1 with inexact real random numbers the sum of whose
  97. ;;squares is less than 1.0.  Thinking of @1 as coordinates in
  98. ;;space of dimension @var{n} = @code{(vector-length @1)}, the
  99. ;;coordinates are uniformly distributed within the unit @var{n}-shere.
  100. ;;The sum of the squares of the numbers is returned.
  101. (define (random:hollow-sphere! vect . args)
  102.   (let ((ms (sqrt (apply random:normal-vector! vect args))))
  103.     (do ((n (- (vector-length vect) 1) (- n 1)))
  104.     ((negative? n))
  105.       (vector-set! vect n (/ (vector-ref vect n) ms)))))
  106.  
  107.  
  108. ;;; For the uniform distribution on the solid sphere, note that in
  109. ;;; this distribution the length r of the vector has cumulative
  110. ;;; distribution r^n; i.e., u=r^n is uniform [0,1], so r can be
  111. ;;; generated as r=u^(1/n).
  112.  
  113. ;;@args vect
  114. ;;@args vect state
  115. ;;Fills @1 with inexact real random numbers the sum of whose
  116. ;;squares is equal to 1.0.  Thinking of @1 as coordinates in space
  117. ;;of dimension n = @code{(vector-length @1)}, the coordinates are
  118. ;;uniformly distributed over the surface of the unit n-shere.
  119. (define (random:solid-sphere! vect . args)
  120.   (apply random:hollow-sphere! vect args)
  121.   (let ((r (expt (random:uniform1 (if (null? args) *random-state* (car args)))
  122.          (/ (vector-length vect)))))
  123.     (do ((n (- (vector-length vect) 1) (- n 1)))
  124.     ((negative? n))
  125.       (vector-set! vect n (* r (vector-ref vect n))))))
  126.